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SUMMARY 

A  two-step  approach  to  finite  element  ordermg  is  introduced.  The  scheme  involves  ordermg  of 
the  finite  elements  first,  based  on  their  adjacency,  followed  by  a  local  numbering  of  the  nodal 
variables.  The  ordermg  of  the  elements  is  performed  by  the  Cuthill-McKee  algorithm.  This 
approach  lakes  into  consideration  the  underlying  structure  of  the  finite  element  mesh,  and  may  be 
regarded  as  a  "natural"  finite  element  ordering  scheme.  The  experimental  results  show  that  this 
two-step  scheme  is  more  efficient  than  the  reverse  Cuthill-McKee  algorithm  applied  directly  to  the 
nodes,  m  terms  of  both  execution  time  and  the  number  of  fill-in  entries,  particularly  when  higher 
order  finite  elements  are  used.  In  addition  to  its  efficiency,  the  two-step  approach  increases 
modularity  and  flexibility  in  finite  element  programs,  and  possesses  potential  application  to  a 
number  of  finite  element  solution  methods. 

INTRODUCTION 

The  use  of  finite  element  methods  typically  involves  solving  a  system  of  equilibrium  equations  : 
K  u  *  R  (1) 

where  u  and  R  are,  respectively,  the  displacement  and  the  loading  vectors.  K  is  the  global 

stiffness  matrix  which  is  often  symmetric,  positive  defimte,  and  populated  with  many  zeros.  In 

the  solution  process,  it  is  important  to  take  into  account  the  structure  of  the  matrix  K.  i.e.  the 

pattern  of  the  zero  and  nonzero  entries  so  that  the  computational  effort  is  minimized. 

The  solution  process  for  solving  Equation  1  can  be  divided  into  four  separable  tasks'*  : 

1.  Ordermg  -  to  find  a  proper  permutation  matrix  P,  such  that  the  symmetric  permuted 
matrix  PKP‘  has  a  desirable  structure. 

2.  Storage  allocation  -  to  determine  the  necessary  information  about  the  structure  of  the 
matrix  factor  of  PKP‘  and  to  set  up  the  storage  scheme. 

3.  Numeric  factorization  -  to  decompose  the  permuted  matrix  PKP'  into  a  triple  matrix 
product  LDL‘. 

4.  Solution  -  to  compute  the  solution  vector  u  successively  by  a  forward  substitution, 
z  •  L''(P‘R),  and  a  backward  substitution,  u  «  PILL’D*' z). 

The  identification  of  these  four  separable  tasks  not  only  encourages  software  modularity,  but  also 
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facilitates  tbe  theoretical  study  of  sparse  matrices.  This  paper  focuses  on  the  first  task,  that  of 
ordering  of  a  system  of  equilibrium  equations  resulting  from  finite  element  discretization. 

When  the  matrix  K  is  decomposed  into  its  matrix  factors  L  and  D.  it  typically  suffers  some 
fill-in;  that  is,  the  filled  matrix,  F  «  L  L‘,  has  nonzeros  in  positions  which  are  zero  m  K.  One 
objective  in  ordering  a  matrix  is  to  minimize  the  number  of  these  fill-in  entries.  By  doing  so, 
the  computing  cost  may  also  be  reduced  since  the  amount  of  computation  depends  on  the  number 
of  nonzeros  m  the  matrix  factor  L. 

A  finite  element  mesh  consists  of  a  collection  of  finite  elements  which  are  connected  at  their 
common  boundaries.  The  discretization  of  the  mesh  and  the  finite  elements  are  selected  so  that 
the  behavior  of  the  physical  structure  under  examination  can  be  properly  modelled.  The  nodal 

variables  are  then  defined  on  the  finite  element  mesh;  their  number  depends  on  the  type  of  finite 

elements  used.  Traditionally,  most  ordering  schemes  have  been  developed  for  ordering  directly  the 
nodal  variables  and,  hence,  the  equilibrium  equations  in  1.  Very  little  effort  has  been  reported  m 
developing  ordering  schemes  which  include  consideration  of  the  underlying  topological  structure  of 
the  finite  element  mesh. 

One  solution  scheme  specifically  designed  for  solvmg  a  system  of  equations  resultmg  from  finite 
element  problems  is  tbe  frontal  solution  method*'*.  This  method  is  often  considered  as  tbe  most 
"natural"  solution  scheme  since  it  operates  directly  on  the  underlying  structure  of  tbe  finite 

element  mesh.  In  this  solution  scheme,  the  finite  elements  are  assembled  and  entered  mto  the 
solution  one  at  a  time.  A  nodal  variable  is  eliminated  as  soon  as  all  the  neighboring  finite 

elements  mcident  to  it  are  available.  The  remainmg  assembled  nodal  variables,  which  are  not  yet 
ready  for  elimination,  are  retained  in  the  "front",  or  "active  front". 

Another  solution  scheme  which  is  similar  to  the  frontal  solution  method  is  the  so-called 
generalized  element  method^°  In  this  method,  the  nodal  variables  of  an  active  front  are  treated 
as  a  superelement  Multiple  fronts  are  allowed,  and  each  front  is  treated  as  a  superelement. 
This  method  has  recently  attracted  attention  from  numerical  analysts^.  The  merit  of  the  frontal 
solution  and  the  generalized  element  methods  is  that  both  can  be  extended  to  include 
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considerations  of  auxiliary  storage  in  a  simple  manner. 

Obviously,  the  efficiency  of  the  frontal  solution  and  the  generalized  element  methods  depends  on 
the  order  in  which  the  finite  elements  are  numbered.  The  ordermg  of  the  nodal  variables  is 
implicitly  imposed  by  the  order  in  which  the  finite  elements  are  processed.  To  the  authors' 
knowledge,  there  have  not  been  any  studies  reported  for  ordering  the  finite  elements  directly.  In 
this  paper,  we  attempt  to  examine  the  feasibility  of  ordering  the  finite  elements. 

It  will  be  shown  that  ordering  the  finite  elements  has  considerable  computational  advantage  over 
ordering  the  nodal  variables,  but  that  ordering  the  finite  elements  alone  does  not  completely  define 
the  permutation  matrix  P.  It  is  therefore  necessary  to  include  into  the  ordering  scheme  a  local 
ordering  strategy  for  numbering  the  nodal  variables  of  each  finite  element  so  that  the  number  of 
fill-in  entries  can  be  minimized.  This  "two-step”  approach,  which  includes  the  ordermg  of  the 
finite  elements  as  well  as  of  the  nodal  variables,  is  the  subject  of  this  paper. 

GRAPH  THEORY  NOTATION  AND  SPARSE  MATRIX  STUDY 
The  use  of  graph-theoretic  approaches  has  found  many  applications  in  sparse  matrix  study. 
Although  few  results  from  graph  theory  are  directly  applicable  to  the  study  of  sparse  matrices, 
the  graph  representation  is,  nonetheless,  a  powerful  tool  to  characterize  the  structure  of  a  sparse 
matrix. 

Symmetric  Matrices  and  Finite  Graphs 

A  finite  undirected  graph  G=(V,£)  consists  of  a  finite  set  V  of  n  elements,  called  nodes,  and 

a  set  E  of  unordered  pairs  of  distinct  nodes  (v.v).  called  edges.  For  any  node  v  in  G,  the  set 

>  j 

of  nodes  adjacent  to  v,  adjtv),  is  defined  as 
adj(v)  =  {  v  (  V  I  (v.v  )  <  E  }. 

i  i 

The  number  of  nodes  in  adj(v),  denoted  by  ladj<v)l,  is  called  the  degree  of  v.  The  deficiency 
of  v,  def(v),  is  the  set  of  distinct  pairs  of  nodes  in  adj<v)  which  are  not  themselves  adjacent.  A 
graph  is  complete  if  every  pair  of  nodes  is  adjacent.  A  subgraph  G'=(V’,E’)  of  G  is  a  graph 
for  which  V’  c  V  and  E'  c  E.  A  complete  subgraph  is  called  a  clique.  A  section  graph 
G(V’)  IS  a  «ibgraph  G*(V’,E(V’))  induced  by  a  node  set  V,  where 
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E(V’)  =  {  (v,v)  «  E  1  v.v  (  V’} 

>  j  >  j 

If  a  set  of  nodes  V"  is  deleted  from  the  graph  G*(V,E),  the  section  graph  G(V\V")  is  obtained 

from  G  by  removing  the  nodes  V"  together  with  their  incident  edges. 

A  (simple)  path  {v,  v.}  is  a  sequence  of  distinct  nodes  and  continuous  edges  leading  from 

V  to  V  such  that  there  are  no  repeating  edges.  A  graph  G  is  said  to  be  connected  if  there  is 
at  least  one  path  between  every  pair  of  distinct  nodes  in  G;  otherwise,  G  is  disconnected.  A 
connected  subgraph  is  called  a  component.  The  distance  d(  v  ,v )  between  two  nodes  v  and  v 

■  J  ‘  J 

in  a  connected  graph  is  the  number  of  edges  in  the  shortest  path  joining  nodes  v  and  v.  The 

eccentricity  e(v)  of  a  node  v  is  then  given  by 
e(v)  »  maa  {  d(v,v )  1  v  «  V  }. 

i  i 

The  diameter  5(C)  of  the  graph  G=(V.E)  is  defined  as 
5(C)  •  max  {  e(v),  v  <  V  }. 

A  node  v  is  said  to  be  a  peripheral  node  if  e(v)  =  5(C). 

For  a  graph  G»(V,E)  with  n  nodes,  an  ordering  (or  labelling)  of  V  is  a  bijection  (a  one-to- 

one  onto  mapping)  ; 

a  :  i  I,  2 .  n  }  <->  V. 

The  ordered  graph  of  G  is  denoted  by  G^  *  (V,E,«).  The  integer,  rangmg  from  1  to  n, 

assigned  to  a  node  by  the  ordering  is  called  the  number  or  label  of  that  node. 

The  notation  of  a  level  structure  is  commonly  used  to  describe  the  properties  of  many 

ordermg  schemes.  A  level  structure  LS  of  a  graph  G  is  an  arrangement  of  the  nodes  of  G  into 

m  levels,  ordered  LS  ,  ....  LS  ,  such  that  nodes  at  a  given  level  LS  are  connected  to  no  nodes 

other  than  LS  ,  LS  and  LS  For  a  node  v  «  V,  there  corresponds  a  level  structure  rooted 

k-l  k  k-1 

at  V.  The  levels  of  the  rooted  level  structure  are  determined  by 

1.  assigning  LS^  =  {v},  and 

2.  for  k  >  1,  assigning  to  LS^  the  set  of  nodes  adjacent  to  nodes  m  LS^  ^  which  have 
not  yet  been  included  in  any  previous  levels. 

Given  an  n  by  n  symmetric  matrix  K.  there  corresponds  to  it  a  finite  graph  G*(V.E).  where  a 
node  V  «  V  denotes  the  i"'  row  (or  column)  of  the  matrix  K.  and  an  edge  (v,v>  t  E 

'  '  '  J 
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symbolizes  an  offdiagonal  nonzero  entry  K  (=  K  ).  An  ordermg  of  the  graph  C  of  matrix  iC  is 
equivalent  to  a  symmetric  permutation  PKP‘.  where  P  is  an  n  by  n  permutation  matru.  The 
unordered  graphs  (in  which  nodes  are  not  labelled)  of  K  and  PKP‘  are  the  same,  but  the  node 
numbers  of  the  associated  ordered  graphs  are  different. 

Finite  Element  Mesh  and  Its  Solution  Graph 

A  finite  element  mesh  is  a  collection  of  finite  elements  m  which  the  adjacent  finite  elements 
are  joined  at  their  common  boundaries  or  vertices.  There  is  a  node  at  each  vertex  of  the  finite 
element  mesh.  For  finite  elements  of  higher  order,  where  higher  order  polynomial  mterpolation 
functions  are  used,  there  may  also  be  nodes  lying  along  the  sides  or  on  the  faces  of  the  fimte 
elements,  and/or  located  internally  within  the  finite  element  itself.  The  nodes  situated  at  the 
interior  of  the  finite  element  are  referred  as  the  interna/  variables,  and  the  nodes  located  along 
the  boundaries  are  called  the  externa!  variables. 

Associated  with  each  node  is  a  set  of  variables  corresponding  to  the  degrees  of  freedom  defined 
at  that  node.  This  set  of  variables,  in  turn,  corresponds  to  a  submatrix  m  the  global  stiffness 
matru.  which  is  typically  full,  Durmg  assembly  and  solution,  this  submatrix  can  conveniently  be 
treated  as  a  smgle  entity.  For  simplicity,  the  set  of  variables  at  a  node  is  referred  as  the  nodal 
variable. 

The  finite  element  mesh  can  be  transformed  directly  into  the  fmite  graph  G  representmg  the 
structure  of  the  global  stiffness  matrix.  We  call  this  graph  a  finite  element  solution  graph. 
or,  simply,  a  solution  graph  (to  distmguish  it  from  the  finite  element  connectivity  graph  to  be 
introduced  below).  The  nodes  of  the  solution  graph  are  the  nodal  variables  defined  on  the  finite 
element  mesh.  The  edges  of  the  solution  graph  are  constructed  by  makmg  the  nodes  of  each 
finite  element  pairwise  adjacent,  smce  a  nonzero  entry  K  in  the  global  stiffness  matru  K  implies 
that  the  i***  and  the  j"*  nodal  variables  share  at  least  one  incident  finite  element.  In  other  words. 


if  V’  denotes  the  set  of  nodal  variables  belonging  to  one  finite  element,  the  section  graph  G(V') 
of  the  solution  graph  G  is  a  clique:  the  union  of  the  cliques  over  all  elements  defines  the  edges 
m  the  solution  graph.  Examples  of  transforming  a  fimte  element  mesh  into  its  associated  solution 


6 


graph  are  shown  m  Figure  1. 

Matrix  Factorizatioo  and  Its  Graph  Theoretic  Model 

The  factorization  of  K  can  be  symbolically  modelled  using  a  graph-theoretic  approach.  This 

approach  is  particularly  helpful  m  understanding  how  the  fill-in  entries  are  created  durmg 

factorization. 

The  most  fundamental  scheme  of  matrix  factorization  is  one  analogous  to  Gaussian  elimination, 
summarized  in  Figure  2.  At  the  i'*"  step  of  the  factorization,  the  i'**  column  of  L  and  the  i'*" 
diagonal  entry  of  D  are  computed,  and  the  matrLX  K*'*  is  condensed  and  modified  to  It  is 

in  the  condensation  that  the  fill-in  entries,  if  any.  are  created. 

Usmg  the  finite  graph  representation  of  a  sparse  matrix,  the  factorization  scheme  can  be 

modelled  graph-theoretically  as  a  node-elimination  process’^  Upon  elimination  of  node  v.  the 

elimmation  graph  G  is  obtamed  from  G-<V.£)  by  : 

1.  deleting  node  v  and  its  incident  edges: 

2.  adding  auxiliary  (fill-m)  edges  such  that  all  adjacent  nodes  of  v  form  a  clique. 

That  IS. 

G  •  (V\v,  E(V\v)  U  def(v)). 

V 

For  an  ordered  graph  G^,  the  elimination  is  a  recursive  process  defined  by 
P(G„)  =■  {  G  '  G„.  G .  G  } 

a  0  «  1  n-l 

where  G  is  called  the  i'*"  elimination  graph  obtamed  bv  eliminating  v  from  G  The  fill-in 
cdees  created  durine  the  elimination  of  v,  denoted  bv  f,  are  the  def(v  )  in  G  These  fill-m 

-  -  i  '  1  11-1 

edges  correspond  to  the  new  nonzero  entries  introduced  mto  the  condensed  matrix  K"*"  at  the  i'*' 
step  of  the  factorization. 

In  a  f'mte  element  soiution  graph,  the  new  clique  formed  by  the  elimination  of  a  node  can  be 
treated  as  a  new  finite  element,  called  a  generalized  element  or  a  superelement"®.  The  node- 
elimination  process  on  a  finite  element  graph  can  thus  be  interpreted  as  a  series  of  transformation 
of  the  fimte  elements  into  the  superelements.  This  elimination  model  also  suggests  an  mteresting 
characteristic  related  to  higher  order  finite  elements.  Internal  nodes  of  a  higher  order  element  are 


j 
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connected  only  to  the  external  nodes  of  that  element.  Furthermore,  all  nodes  m  a  finite  element 
form  a  clique.  Hence,  if  an  internal  node  is  eliminated  before  any  of  the  external  nodes  of  the 
same  element,  there  is  no  fill-in  smce  the  deficiency  of  the  mternal  node  is  null. 

After  elimination,  the  set  of  fill-m  edges  is  given  by 
n-1 

<I><G  )  =  U  f. 

"  1=1  ‘ 

The  filled  graph  G^,  which  represents  the  structure  of  the  filled  matrix,  F  =  L  *  L',  is  defmed 
by 

G  =  (V,  E  U  4>(G  )). 

F  a 

The  node-elimination  model  was  formally  mtroduced  by  Rose’*.  Since  then,  efforts  have  been 
made  to  develop  efficient  algorithm*;  implementmg  the  model’^'  ".  These  algorithms  have  been 
encoded  mto  sparse  matrix  programs  such  as  YMSL*  and  Sparspak’.  The  major  application  of 
these  algorithms  is  in  the  second  task  presented  in  the  Introduction,  namely  to  set  up  the  data 
structure  for  storing  the  numeric  entries  of  the  matrix  factors. 

A  node-addition  model  which  simulates  the  Cholesky  factorization  has  recently  been  developed 
by  the  authors'^.  For  this  graph-theoretic  model,  nodes  are  added  onto  the  filled  graph  G^,  rather 
than  being  elimmated  from  the  original  graph  G^.  The  structure  of  the  matrix  factor  L  is 

constructed  one  row  at  a  time.  When  computmg  the  entries  of  a  particular  row  m  L.  the  model 
does  not  require  any  a  priori  information  for  the  rows  beyond  the  current  row.  Therefore,  this 
model  has  the  flexibility  that  the  tasks  of  labelling  a  nodal  variable,  determinmg  its  row  structure 

in  L.  and  computing  the  numeric  entries  of  the  row  can  all  be  performed  simultaneously. 

The  execution  limes  for  the  numeric  factorization  as  well  as  the  symbolic  factorization  of  a 
given  matrix  fC  depend  on  the  number  of  nonzeros  in  the  matrix  factor  L.  Therefore,  it  is 

worthwhile  to  develop  efficient  algorithms  to  order  the  matru  K  so  that  the  number  of  fill-m 
entries,  or  equivalently  the  number  of  nonzeros  in  L,  is  minimized.  In  the  next  section,  two  well 
known  ordering  algorithms,  namely  the  Cuthill-McKee  algorithm  and  the  reverse  Cuthill-McKee 
algorithm,  are  discussed. 
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The  Cuthill-McKee  and  the  Reverse  Cuthill'McKee  Ordering  Algorithms 

Let  K  be  an  n  by  n  symmetric  matru.  For  the  i"'  row  of  K.  defme 

(7  (K)  =  min  1  j  1  K  ^  0  } 

‘  >J 

and 

*  1  -  O’  (K). 

1  I 

The  number  <r  (K)  denotes  the  column  subscript  of  the  first  nonzero  entry  in  the  1“*  row  of 

matrix  K.  The  number  j3  (K.)  is  usually  referred  as  the  local  bandwidth  of  the  row  of 

matrix  K,  and  is  equal  to  the  number  of  off-diagonal  entries  between  the  first  nonzero  of  the 

row  and  the  mam  diagonal.  The  bandwidth  and  the  profile  of  matrLX  K  are  defined  as 
b(K)  =  max  {  ^(K),  i  *  1,  ...n  } 

and 

n 

/)<K)  =  I  /ff(K) 

1=1  ’ 

respectively.  The  nonzero  entries  of  the  filled  matrix  F  of  K  are  confmed  withm  the  local  band 
of  each  row.  Many  ordering  schemes  have  been  developed  to  mimmize  the  bandwidth  or  the 
profile  of  K. 

One  ordering  scheme  commonly  used  with  finite  element  programs  to  minimize  the  bandwidth  is 
the  Cuthill-McKce  algorithm'*.  Basically,  the  algorithm  is  a  breadth-first  technique  in  labellmg 
the  nodes  of  a  graph.  Once  a  finite  element  mesh  is  tranformed  into  its  finite  element  solution 
graph,  the  algorithm  numbers  the  nodes  in  the  following  way. 

1.  Determine  a  starting  node  and  number  it  v^. 

2.  For  nodes  v,  i=l,...,n.  find  all  unlabelled  neichborinc  nodes  of  the  node  v  and 
number  them  sequentially  in  mcreasmg  order  of  degree. 

This  ordermg  algorithm  is  essentially  the  same  as  that  for  generating  a  level  structure  rooted  at 
the  startmg  node  vy 

In  selecting  a  starting  node,  a  peripheral  node  is  preferred'^.  The  idea  of  choosing  a  peripheral 
node  as  a  starting  node  is  to  generate  as  many  levels  as  possible  in  the  level  structure,  smce  an 
mcrease  in  the  number  of  levels  tends  to  decrease  the  profile  of  the  corresponding  matrix. 
Unfortunately,  existing  algorithms  for  finding  a  peripheral  node  are  not  computationally  feasible. 
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A  compromise  that  has  been  suggested  is  to  choose  a  starting  node  H'lth  high  eccentricity.  This 
node  is  generally  referred  as  a  pseudo- peripheral  node.  An  efficient  algorithm  to  find  a 
pseudo-peripheral  node  has  been  suggested  and  implemented  by  George  and  Liu'°.  For  most 
practical  cases,  the  execution  time  for  finding  a  pseudo-peripheral  node  is  no  greater  than 
0<,'£|),  i.e.,  the  order  of  the  number  of  edges  m  the  graph  G=(V,E). 

George  eUal.*-  have  shown  that  if  linear  insertion  is  used  for  sortmg  the  time  complexity  for 

the  second  step  of  the  Cuthill-McKee  algorithm  requires  ai  most 
(  4 1  £  1  +  2cm  I  £  I  )  operations, 

where  l£!  is  the  number  of  edges  m  the  graph  G,  m  is  the  maximum  degree  of  any  node,  and  c 
is  some  constant. 

It  has  been  discovered  that  significant  improvements  in  minimizing  the  profile  can  be  achieved 
by  simply  reversing  the  node  ordering  obtained  from  the  Cuthill-McKee  algorithm*.  The  resultmg 
algorithm  is  the  well  known  reverse  Cuthill-McKee  (RCM)  algorithm,  which  can  be  summarized 
as  follows  : 

1.  number  the  nodes  by  the  Cuthill-McKee  algorithm:  and 

2.  reorder  the  nodes  by  reversing  the  node  numbering  obtamed  m  Step  1. 

To  reverse  the  ordermg  of  n  nodes  requires  only  n  operations.  Therefore,  the  overall  complexity 
of  the  RCM  algorithm  remains  bounded  by  0<mi£!). 

It  has  been  proved  by  Liu  and  Sherman  that  the  RCM  algorithm  is  never  mferior  to  the 
Cuthill-McKee  algorithm'^.  In  its  application  to  finite  element  ordering,  they  also  found  that  the 
RCM  algorithm  is  particularly  superior  when  finite  elements  of  higher  orders  are  used. 

A  listmg  of  computer  subroutines  for  the  RCM  algorithm  can  be  found  in  reference  12.  A 
more  detailed  description  of  the  algorithm  can  also  be  found  in  that  reference. 
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A  TWO-STEP  APPROACH  TO  FINITE  ELEMENT  ORDERING 

In  Hus  section,  we  present  a  "two-step”  approach  to  finite  element  ordermg.  The  ordermg 
process  is  divided  into  two  separate  tasks.  The  first  task  orders  the  finite  elements  in  the  finite 
element  mesh.  The  second  task  then  labels  the  nodal  variables. 

Representation  of  Finite  Element  Connectivity 

The  method  to  be  described  requires  an  explicit  representation  of  the  connectivity  of  the  finite 
elements  in  the  mesh.  The  connectivity  of  the  fmite  elements  can  be  topologically  represented  as 
a  graph;  we  call  this  graph  a  finite  element  connectivity  graph,  or.  simply,  a  connectivity 
graph.  The  nodes  m  the  connectivity  graph  correspond  to  the  finite  elements  m  the  mesh.  The 
edges  of  the  connectivity  graph  are  used  to  describe  the  interconnections  between  the  fimte 

elements.  One  possible  way  to  define  the  interconnections  between  the  fimte  elements  is  to  say 
that  two  finite  elements  are  adjacent  if  they  share  a  common  node  in  the  mesh.  This  defmition. 
however,  may  lead  to  a  huge  number  of  edges  m  the  connectivity  graph.  Since  the  efficiency  of 
most  existing  ordering  algorithms  is  a  function  of  the  number  of  edges  in  the  graph,  this 
representation  of  fimte  element  connectivity  may  suffer  a  significant  drawback  in  terms  of  the 
execution  tune  required  for  ordering  the  finite  elements.  Here,  we  examine  an  alternative 
definition  m  which  the  number  of  edges  in  the  finite  element  connectivity  graph  can  be  vastly 
reduced. 

In  finite  element  methods,  a  continuum  is  subdivided  by  imaginary  lines  or  surfaces  into  a 
number  of  finite  elements.  These  elements  are  assumed  to  be  mterconnected  at  the  vertices 
situated  on  their  boundaries.  Therefore,  it  is  simple  to  assert  that  the  finite  elements  are 
topologically  interconnected  by  their  boundaries.  A  fimte  element  of  n  dimensions,  where 

n  =  1.2.3.  posseses  boundaries  of  (n-l>  dimensions.  For  instance,  the  boundaries  of  a  3- 

dimensional  (volumetric)  finite  element  are  the  2-dimensional  boundary-faces,  so  that  m  a  3- 
dimensional  contmuum,  the  volumetric  finite  elements  are  mterconnected  by  their  boundary-faces. 
In  the  same  fashion.  2-dimensional  (planar)  finite  elements  are  bounded  and  interconnected  by 
their  1 -dimensional  boundary-lines,  and  1 -dimensional  (linear*  fmite  elements  are  connected  to 
their  adjacent  elements  through  the  0-dimensional  boundary-nodes.  Hence,  m  a  finite  element 


mesh,  two  fimte  elements  are  said  to  be  adjacent  if  they  are  connected  at  their  boundaries,  rather 
than  at  their  vertices.  Usmg  this  defimtion,  the  nodes  in  the  fimte  element  connectivity  graph  are 
the  finite  elements  in  the  fimte  element  mesh,  while  the  edges  connecting  the  nodes  in  the 
connectivity  graph  correspond  to  the  imaginary  boundaries  of  the  fimte  elements  separating  the 
continuum.  Examples  of  finite  element  meshes  and  their  associated  connectivity  graphs  are  shown 
in  Figure  3.  For  a  two-dimensional  (planar)  continuum,  it  is  interesting  to  note  that  the  fimte 
element  connectivity  graph  is  analogous  to  the  dual  of  a  planar  graph^,  with  the  exterior  node  in 
the  dual  graph  omitted.  In  general,  the  element-element  connectivity  relationship  is  topologically 
equivalent  to  the  definition  of  a  dual  of  an  n-dimensional  complex'^- 

In  some  cases,  the  connectivity  of  finite  elements  cannot  be  completely  represented  usmg  the 
definition  given  above.  As  shown  in  Figure  4(a).  n-dimensional  finite  elements  may  not 
necessarily  be  connected  through  all  their  (n-l>-dimensional  boundaries.  Another  example  may  be 
that  the  adjacent  fimte  elements  do  not  have  the  same  geometric  dimensions,  as  illustrated  m 
Figure  4(b).  Nevertheless,  the  transformation  process  described  is  still  valid  and  applicable  to 
these  cases;  the  resulting  finite  element  connectivity  graph  may  at  worst  become  a  disconnected 
graph.  Each  connected  component  in  the  connectivity  graph  can  be  labelled  mdependently  by  the 
method  presented.  The  nodes  at  the  interface  between  two  connected  components  must  then  be 
numbered  higher  than  all  other  nodes  in  the  two  components.  The  e.xamples  just  described  are  not 
very  common  in  practice,  and  will  not  be  pursued  further. 

There  are  several  major  advantages  in  defming  the  connectivity  of  the  finite  elements  by  means 
of  their  boundaries,  instead  of  their  nodes.  First,  the  number  of  edges  in  the  connectivity  graph 
is  minimized.  Furthermore,  this  defimtion  completely  disregards  the  nodal  variables,  or  vertices, 
in  the  finite  element  mesh.  Therefore,  the  ordering  of  finite  elements  can  be  performed  before 
the  types  of  finite  elements,  and  thus  the  number  of  nodal  variables,  are  chosen. 
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Ordering  the  Finite  Elements 

Once  the  adjacency  structure  of  a  connectivity  graph  has  been  established,  the  Cuthill-McK.ec 
algorithm  can  be  applied  to  number  the  nodes  of  the  connectivity  graph,  which  correspond  to  the 
fmite  elements  m  the  mesh.  An  example  of  a  S  by  S  regular  mesh  with  triangular  fmite 

elements  is  shown  in  Figure  5(a).  The  mesh  is  first  transformed  mlo  its  connectivity  graph 

shown  in  Figure  5  o).  The  nodes  m  the  connectivity  graph  are  ordered  by  the  Cuthill-McKec 

algorithm.  The  resultmg  ordering  of  the  fimte  elements  is  given  in  Figure  5(c). 

As  noted,  earlier,  the  time  comple.xity  of  the  Cuthlll-McKee  algorithm  is  a  function  of  the 

number  of  edges  m  a  graph.  For  a  given  fimte  element  mesh,  the  number  of  edges  m  the 

associated  connectivity  graph  is  considerably  less  than  the  number  of  edges  m  its  finite  element 
solution  graph  defined  previously.  When  higher  order  finite  elements  are  used,  the  difference  is 
even  more  dramatic,  smce  the  number  of  edges  m  a  finite  element  solution  graph  increases 

combinatonally  with  the  number  of  nodal  variables  per  fimte  element.  On  the  other  hand,  the 
number  of  edges  in  the  finite  element  connectivity  graph  remams  constant  for  a  given 

discretization  of  the  fmite  element  mesh. 

Let  the  fimte  element  mesh  consist  of  ne  finite  elements,  and  let  nb  denote  the  number  of 

1 

boundaries  of  finite  element  i.  The  overall  time  complexity  of  the  Cuthill-McKee  algorithm, 
following  from  the  previous  discussion,  is  0<m!E!).  where  m  is  the  largest  value  of  nb . 

I 

1  -  1.  ....  ne,  and  lEt  is  the  total  number  of  internal  boundaries  interconnecting  the  fimte 

elements. 

Ordering  the  Nodal  Variables 

A  good  ordering  of  the  fimte  elements  does  not  automatically  guarantee  a  good  node  ordering, 
m  the  sense  that  the  number  of  fill-m  entries  is  minimized.  The  ordermg  of  the  finite  elements 
does  not  complete  the  ordering  of  K.  as  the  nodal  variables  associated  with  the  element  do  not 
acquire  labels  by  this  process.  Therefore,  one  must  also  consider  the  ordermg  of  the  nodal 
variables.  The  proposed  strategy  is  to  label  the  nodal  variables  of  the  finite  elements  following 
the  order  in  which  the  finite  elements  are  numbered.  For  each  fmite  element,  the  nodal  variables 
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are  ordered  according  to  their  valencies,  where  the  valency  of  a  node  is  the  number  of  finite 
elements  incident  at  that  node.  The  nodal  valency  is  an  mdicator  of  the  degree  of  a  node. 

The  local  ordering  scheme  to  label  the  nodal  variables  is  summarized  as  follows  ; 

1.  Determine  the  nodal  valency  for  each  node  in  the  finite  element  mesh. 

2.  (Mam  loop)  Enter  the  finite  elements  one  at  a  time  following  the  order  m  which  the 
finite  elements  were  previously  numbered  by  applying  the  (3uthill-McKee  algorithm  to 
the  connectivity  graph.  For  each  finite  element,  find  all  unlabelled  nodes  connected 
to  It  and  number  them  sequentially  in  descreasing  order  of  their  valencies. 

3.  (Reverse  ordermg)  Reverse  the  ordering  of  the  nodal  variables  obtamed  m  Step  2. 

In  Step  2  of  the  local  ordermg  scheme,  a  nodal  variable  with  mmimum  valency  among  the 
unlabelled  nodal  variables  in  a  fimte  element  is  numbered  last  The  node  ordermg  obtamed  is 
then  reversed  m  Step  3.  The  motivation  behind  these  two  steps  is  that  the  nodal  variables  with 
the  lowest  valency  in  a  finite  element  will  be  eliminated  first  This  strategy  of  minimum  degree 
(valency)  ordermg  has  been  employed  in  many  popular  node  ordering  schemes,  such  as  the 
mmimum  degree  ordermg  algorithm'*  and  the  Cuthill-McK.ee  algorithm'*,  although  the  strategy  is 
used  in  a  different  manner  in  different  schemes. 

It  has  been  mentioned  that  by  reversing  the  node  ordering  obtained  from  the  Cuthill-McKee 

algorithm,  the  result  can  be  improved  significantly.  This  property  is  also  true  in  the  local 

ordering  scheme  proposed.  First,  reversing  the  numbermg  of  the  nodal  variables  which  are 
ordered  in  decreasing  order  of  their  valencies  ensures  that  internal  nodal  variables  are  numbered 
before  the  external  variables  of  the  same  element,  so  that  the  internal  nodal  variables  are 
eliminated  before  the  external  nodal  variables  of  the  same  finite  element.  Hence,  there  will  be  no 
fill-in  created  when  elimmatmg  the  internal  nodal  variables.  This  ordering  strategy  also  has  the 
property  that  the  nodes  situated  along  a  boundary  are  numbered  and  eliminated  before  the  comer 
nodes  of  the  same  boundary.  This  is  because  the  corner  nodes  are  m  general  connected  to  more 

finite  elements  than  the  nodes  located  along  a  boundary.  The  nodal  numbering  produced  by  the 

proposed  scheme  for  a  simple  example  of  a  5  by  5  regular  mesh  consisting  of  10-node  triangular 
finite  elements  is  shown  m  Figure  6.  The  two  properties  just  described  are  clearly  demonstrated 
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by  this  example. 

For  step  1,  given  tie  element-node  incidence  table  —  a  listing  of  the  nodes  incident  on  the 

fimte  elements  ~  the  determination  of  the  nodal  valencies  can  be  done  in  exactly 
ne 

ITI  »  Z  nv  operations, 

1*1  ‘ 

where  nv  is  the  number  of  nodal  variables  of  the  i*”  finite  element,  ne  is  the  number  of  finite 

i 

elements  in  the  mesh,  and  1 T I  is  the  size  of  element-node  mcidence  table. 

To  implement  step  2,  a  Imear  insertion  may  be  used  to  sort  the  unlabelled  variables  of  each 

fimte  element.  For  some  constant  c,  the  sortmg  of  nv  elements  by  linear  insertion  requires 

c(nv-)  operations’.  Thus,  the  time  complexity  of  Step  2  of  the  algorithm  requires  al  worst 
ne  ne 

c<  Z  nv*)  <  c<  Z  nv )  nv  =  c'TInv  operations, 

.  ,  I  ~  ,  i  max  max 

1=1  1*1 

where  nv  *  max<nv ),  i=l,...,ne.  The  number  of  operations  required  to  reverse  the  node 

mail  i 

ordering  equals  to  the  number  of  nodal  variables,  nn,  in  the  finite  element  mesh.  Hence,  the 
time  complexity  for  numbering  the  nodes,  given  the  ordering  of  the  fmite  elements,  is  bounded  by 
(0(  ITlnv  )  *  nn). 

max 


EXPERIMENTAL  RESULTS 

In  this  section,  we  report  some  experimental  results  comparing  the  reverse  Cuthill-McKee 
algorithm  and  the  two-step  ordermg  algorithm.  Four  fimte  element  models  have  been  selected  as 
examples.  They  are  :  (1)  a  5  by  S  regular  mesh  with  triangular  fmite  elements:  (2)  an  L- 
shaped  mesh  with  triangular  finite  elements;  (3)  a  cross-shaped  mesh  with  recUngular  fimte 
elements:  and  (4)  a  2  by  2  by  2  model  consistmg  of  rectangular  block  fimte  elements.  These 
models  are  shown  in  Figure  7.  For  each  model.  Imear.  quadratic  and  cubic  fmite  elements  have 
been  used.  There  are  altogether  twelve  test  cases. 

The  computer  subroutines  for  the  RCM  algorithm  have  been  adopted  from  the  Sparspak  package 
developed  by  George  and  Liu  at  the  University  of  Waterloo’  and  listed  in  Reference  12.  These 
subroutines  mclude  findmg  a  pseudo-peripheral  node  of  a  graph  and  ordering  the  nodes  by  the 
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RCM  algorithm.  The  input  information  to  these  subroutines  is  the  number  of  nodal  variables  and 
an  adjacency  structure  pair  representing  the  node  connectivity  of  the  solution  graph.  A  node 
adjacency  structure  is  illustrated  in  Figure  8U). 

For  the  two-step  ordering  scheme,  the  same  set  of  computer  subroutmes  from  Sparspak  has  been 
employed  to  number  the  finite  elements  by  the  Cuthill-McKee  algorithm,  except  that  the  step  for 
reverse  ordering  has  been  omitted.  In  this  case,  the  input  data  for  the  adjacency  structure  pair  is 
the  node  connectivity  of  the  fmite  element  connectivity  graph.  An  example  of  this  data  structure 
IS  given  m  Figure  8<b).  Once  the  finite  elements  have  been  ordered,  the  adjacency  structure  can 
be  discarded.  The  subroutme  for  performmg  the  second  step,  that  of  ordering  the  nodal  variables. 
IS  listed  in  Appendix  1.  For  this  subroutme,  the  element-node  incidence  table  is  assumed  to  be 
available. 

The  test  cases  have  been  run  on  a  DECsystem  20  computer  at  Carnegie- Mellon  University.  For 
each  test  case,  the  execution  times  for  the  two-step  ordering  scheme  and  the  RCM  algorithm  are 
reported  in  Table  1.  To  mmimize  timing  errors  due  to  a  multi-programmed  operating  system 
environment,  the  test  cases  have  been  run  when  the  computer  was  lightly  loaded.  As  the  results 
mdicate,  the  two-step  ordermg  scheme  requires  much  less  execution  time  than  the  RCM  node 
ordering  algorithm,  except  for  the  cases  when  linear  triangular  finite  elements  are  employed.  For 
those  cases,  the  execution  times  for  ordering  the  nodes  by  the  RCM  algorithm  and  for  ordering 
the  fmite  elements  by  the  Cuthill-McKee  algorithm  are  approximately  the  same.  The  deficiency  of 
the  two-step  ordering  scheme  is  due  to  the  second  step  of  labelling  the  nodal  variables.  For 
cases  when  higher  order  fimte  elements  are  used,  very  large  amounts  of  savings  in  execution 
times  can  be  achieved. 

For  each  lest  case,  the  structures  of  the  matrix  factor  L  resulting  from  the  RCM  and  the  two- 
step  ordermg  algorithms  are  summarized  m  Table  2.  m  terms  of  the  profile  of  the  stiffness 
matrix,  the  number  of  fill-in  entries  and  the  size  of  the  matrix  factor.  For  almost  all  cases,  the 
profiles  resulting  from  the  two-step  ordering  scheme  are  slightly  larger  than  those  obtained  from 
the  RCM  algorithm.  On  the  other  hand,  the  two-step  ordering  scheme,  in  most  cases,  leads  to 
less  fill-in  than  the  RCM  algorithm.  This  result  is  particularly  true  for  finite  element  meshes 


16 


where  higher  order  finite  elements  are  used. 

DISCISSION 

In  this  paper,  a  ''two-step"  finite  element  ordering  scheme  has  'oeen  introduced.  In  addition  to 
the  efficiency  achieved,  the  scheme  is  also  highly  adapta'r'ie  to  various  solution  methods  commonly 
used  in  finite  element  analysts. 

The  major  characteristic  of  the  two-step  ordering  scheme  is  that  ordering  the  finite  elements 

and  ordering  the  nodal  variables  are  separated  into  two  independent  tasks.  T'ne  two  tasks  can  be 
impiemenied  as  two  separate  modules.  This  property  of  modularity  provides  fle\ibiii!y  in  tne 
software  design  of  finite  element  programs.  For  instance,  one  can  choose  to  number  the  finite 

elements  once  the  discretization  of  the  finite  element  mesh  is  established,  e'cen  without  the 
knowledge  about  the  t.vpes  of  finite  elements  to  be  used.  in  view  of  current  soft'^'are 
developments  in  generating  finite  element  meshes  using  graphic  preprocesssors.  t'ne  task  of 

numbering  the  finite  elements  can  easily  be  incorporated  as  an  additional  module  in  the  mesh 
generator  with  very  little  extra  cast. 

The  authors  recognize  that  the  two-step  ordering  scheme  presented  requires  l«'o  sets  of  input 
data:  one  to  descri'oe  the  adjacency  of  the  fimle  elements  and  the  second  for  the  elenicnt-ncde 

incidence  table.  These  two  sets  of  data  can  both  be  generated  by  a  mesh  generator. 
Conventionally,  howes’er.  only  the  element-node  incidence  table  is  created.  In  .\ppendi.\  11.  we 
present  a  computer  program  to  determine  the  element-element  adjacency  structure  pair  using  the 
element-node  incidence  table  as  input.  This  computer  program  ser\es  primarily  for  illustration 
purposes,  and  is  not  meant  to  be  efficient.  The  authors  emphasize  that  the  finite  element 
adjacency  structure  should  be  generated  by  t'ne  mesh  generator,  particularly  w'nen  a  graphical 
preprocessor  is  used. 

Since  the  birth  of  the  finite  element  methods,  efforts  have  been  made  to  deieiop  finite  elements 
using  pol^momial  interpolation  functions  of  higher  orders.  The  e.'.perimenlal  results  favor  strongly 
the  use  of  the  two-step  ordering  scheme  w-uh  higher  order  finite  elements.  In  many  finite 
eieme.nt  problems,  the  finite  elements  are  progressively  modified  by  using  interpolation  functions 
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of  increasingly  higher  order  so  as  lo  improve  the  accuracy  of  the  solution.  At  each  modification, 
the  structure  of  the  global  stiffness  matrix  changes,  but  the  number  and  the  distribution  of  the 
finite  elements  in  the  mesh  often  remain  unchanged.  Examples  of  these  problems  are  quality 
control  in  finite  element  analysis  and  nonlinear  structural  analysis*'.  For  this  type  of  problem, 
the  ordering  of  finite  elements  needs  only  be  determined  once.  The  result  can  be  used  to  reorder 
the  nodal  variables  upon  each  modification  made  to  the  finite  elements.  Moreover,  the  step  m 
ordering  the  nodal  variables  by  the  two-step  ordering  scheme  can  be  performed  much  faster  than 
a  complete  re-numbenng  by  the  node  ordering  schemes. 

In  large  scale  structural  analysis,  substructurmg  is  often  used  to  divide  the  structure  into  two 
or  more  smaller  components,  called  substructures,  which  are  interconnected  at  their  boundaries. 
Topologically,  the  substructures  can  be  treated  m  the  same  manner  as  finite  elements.  Therefore, 
the  strategy  for  ordering  the  finite  elements  may  well  be  applied  to  the  ordering  of  the 
substructures.  In  substructurmg  analysis,  the  external  nodal  variables  located  along  the  boundaries 
of  the  substructures  are  ordered  last.  As  discussed  previously,  this  characteristic  is  also  embedded 
in  the  node  ordering  strategy  of  the  two-step  scheme.  Hence,  there  is  no  reason  that  the  two-step 
approach  caimot  be  extended  to  an  ordering  scheme  for  substructuring  methods. 

One  of  the  characteristics  of  the  frontal  solution  method  is  that  the  assembling  of  the  finite 
elements  and  the  solution  process  are  performed  simultaneously.  For  this  solution  method,  it  is 
the  numbering  of  the  fimte  elements  that  really  matters.  However,  a  proper  ordering  of  the  nodal 
variables  for  each  finite  element  can  further  reduce  the  number  of  fill-in  entries  and,  thus,  the 

computation  cost  of  the  solution  process.  In  the  two-step  scheme  presented,  the  strategy  is  to 

label  the  nodal  variables  following  the  ordering  of  the  finite  elements;  then  the  fmal  node 
ordering  is  reversed.  As  a  result,  the  order  in  w-hich  the  nodal  variables  are  lo  be  eliminated 

ought  lo  follow  the  reverse  order  in  which  the  finite  elements  are  numbered.  Therefore,  by 

entering  the  finite  elements  in  the  reverse  order  of  their  numberings,  the  assembling  and 
factorization  tasks  can  also  be  performed  simultaneously.  However,  unlike  the  frontal  solution 
method,  where  the  nodal  variables  are  arranged  in  an  arbitrary  order,  the  nodal  variables  are  pre¬ 
ordered  by  the  two-step  ordering  scheme  lo  reduce  the  number  of  fill-in  entries.  This  pre- 


ordering  of  the  nodal  variables  has  the  advantage  that  the  data  struoture  for  the  matrix  factors 
can  be  set  up  independently.  Throughout  the  entire  ordering  scheme,  the  structure  of  the  global 
stiffness  matri.v  need  not  e.xtst.  With  the  nodal  variables  pre-ordered.  one  can  proceed  directly  to 
construct  the  data  structure  for  the  matrix  factors  using  e.xisting  ssiiibolic  factorization  algorithms 
Consequently,  while  this  two-step  ordering  scheme  can  improve  the  software  modularity  for 
finite  element  programs,  many  characteristics  of  the  ’‘natural"  frontal  solution  method  can  still  be 
retained  in  the  solution  process. 

The  authors  do  not  claim  that  the  tuo-step  ordering  scheme  proposed  in  this  paper  is  optimal 

in  minimizing  the  number  of  fill-in  entries.  In  fact,  it  has  recently  been  pro\'ed  that  the  problem 

of  computing  the  minimum  fill-in  is  .’^P-complete’’'.  This  ordering  scheme  is  recommended 

because  of  its  efficiency,  modularity  and  fle.xibility.  and  Us  potential  application  to  various  other 
finite  element  problems.  Most  of  all,  the  two-step  ordering  scheme  includes  the  consideration  o! 
the  underlying  topological  structure  of  the  finite  element  mesh,  and  may  therefore  be  regarded  as 
a  "natural"  finite  element  ordering  scheme. 
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Figure  3;  Finite  Element  .Mi>dcl  and  Its  Aseoicaicd 
Finite  Element  Conncctititv  Graph 
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(c)  The  ordcrin"  of  finiic  elements 


Figure  5;  Ordering  the  Finiie  Elcnicnls  bv 
the  Cut  hill  -  McKee  Algorithm 
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(i)  1  inear  element 


(ii)  quadratic  element 


(iii )  cubic  element 


(a)  5  bv  5  regular  mesh  with  iriancular  finite  elements 


(i)  linear  element 


(ii)  quadratic  element 


(iii)  cubic  element 


(b'  L-shpacd  mesh  m  iih  triatiyular  finiit  elements 


.  Fijurc  7;  Finite  Element  Models 
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noce  r.u~ter 


(i)  Finite  element 
solution  graph 
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(iii)  Node-node  adjacency  structure  pair  (IA,JA) 
(a)  Adjacency  structure  pair  of  finite  element  solution  graph 
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(ii)  Element  connection  table 
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(iii)  Element-element  adjacency  structure  pair  (IA,JA) 

(b)  Adjacency  structure  pair  of  finite  element  connectivity  graph 


Figure  8:  Examples  of  Adjacency  Structure  Pair 
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Appendix  1.  Fortran  Computer  Program  for  Labelling  the  Nodal  Variables 

Given  the  ordering  of  the  finite  elements,  the  subroutine  listed  below  labels  the  nodal  vanabl 
in  the  finite  element  mesh. 


C/  »»»»»»»»«**»*»*«»»»»»»»»***»»»*»**»»»*»«»«»»»»»»»»**«»*»»**»»«»»»*»«:»»« 
C/»  • 

C/«  Subroutine  NODORD  :  * 

C/»  Given  t.he  ordering  of  the  finite  elements,  ELPK ,  » 

C/*  this  subroutine  N'CDORD  labels  the  nodal  variables  * 

C/*  as  follows  * 

C/*  1.  Determine  t.he  valency  of  each,  node,  given  the  « 

C/*  element-node  incidence  table  INC  * 

C/*  2.  Label  the  nodes  of  each  finite  ele.ment,  in  decreasing  * 

C/»  order  of  their  valencies  (linear  insertion  is  used)  » 

C/»  3.  Reverse  the  node  ordering  obtained  from  step  2.  » 

c/«  • 

C/*  Input  Variables  :  • 

C/»  INC  -  node-element  INCidence  table  » 

C/»  ELPM  -  ELement  PerKutation  vector  » 

C/«  NV  -  Number  of  Vertices  in  the  finite  element  mesh  » 

C/«  NVE  -  Number  cf  Vertices  per  finite  Element  * 

C/«  NE  -  Number  of  finite  Elements  in  the  mesh  * 

C/«  * 

C/*  Output  Variables  :  « 

C/*  PERK  -  node  PERKutation  vector  * 

C/»  » 

C/*  Woricing  Variables  ;  » 

C/*  INVP  -  iNVrse  Permutation  vector  » 

C/»  NDVL  -  NoDe  Valency  vector  » 

C/»  » 

SUERC’lTINE  NCD0RD(£L?K,  INC,  PER.K,  INVP,  NDVL,  NV,  NVE,  NE) 
IMPLICIT  INTEGER  (A-2) 

DIMENSION  INC(KE,NVE),  ?ER.«(NV),  IKV?(NV1,  EL?.K(NE),  NDVL(NV) 

C/« 

C/*  Initialization 

c/* 

DO  10  K  =  1,NV 
?ERM(K)  =  0 
INV?(K)  =  0 
NDVL(K)  =  0 
10  CONTINUE 

0/ * 

C/»  Determine  node  valency 

C/* 

DO  20  I  =  1,NE 
DO  20  K  =  1,NVE 
NODE  =  INC(:,K) 

IF(NCDE.EQ.O)  GOTO  20 
NDVL (NODE)  =  NDVL (NODE)  +  1 
20  CONTINUE 

C/« 

C/*  Enter  the  elements  sequentially  following  their  orderings 


NXTNCD  =  0 


DO  150  I  =  1,NE 
SLEM  =  ELPM(I) 

START  =  NXTNOD  +  1 
DC  110  K  =  1,NVE 

NODE  =  INC(EEEM,  K) 

IF(INVP(NODE) .NE.O.OR.NODE.EQ.O)  GOTO  110 
NXTNOD  =  NXTNOD  -r  1 
PERM (NXTNOD)  =  NODE 
INVP(NODE)  =  NXTNOD 
110  CONTINUE 

0/» 

C/»  Local  node  ordering  :  nodes  are  sorted  in  decreasing 
C/*  order  of  their  valencies  using  linear  insertion 

c/» 

IF (START  .GE.  NXTNOD)  GOTO  15C 
K  =  START 
120  L  =  K 

K  =  K  +  1 
PERMK  =  PERM(K) 

NDVLK  =  NDVL( PERMK) 

130  IF(L  .LT.  START)  GOTO  140 

PERML  =  PERM(L) 

NDVLL  =  NDVL( PERML) 

IF(NDVLL  .GE.  NDVLK)  GOTO  140 
PERM(L+1)  =  PERML 
L  =  L  *■  1 
GOTO  130 


140 

PERM(L+1) 

=  PERMK 

IF(K  .LT. 

NXTNOD) 

GOTO  120 

150 

CONTINUE 

C/« 

c/« 

Now,  reverse 

the  node 

ordering 

c/« 

KMID  =  NV/2 

J  =  NV 


DO  160  K  =  l.KMID 
I  =  PERK(K) 
PERM(K)  =  PERM(J) 
PERM(J)  =  I 
J  =  J  -  1 
160  CONTINUE 
.RETURN 
END 
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Appendix  II.  Fortran  Computer  Program  for  Generating  the  Adjacency  Structure  Pair  of 
Element -Element  Connectivity 

Given  the  element-node  incidence  table,  the  following  set  of  subroutines  determines  the 
adjacency  structure  pair  of  element-element  connectivity  of  an  arbitrary  finite  element  mesh. 
These  subroutines  are  included  only  for  illustrating  how  the  element-element  adjacency  may  be 
generated  if  the  mesh  generator  produces  only  the  element-node  incidence  table.  The  subroulmes 
are  not  efficient  m  storage,  in  that  space  is  provided  for  the  entire  ssmimetric  node-node 
adjacency  table,  denoted  by  NEL.  Since  this  table  is  very  sparse,  the  function  H.^SK  could  be 
replaced  by  one  that  exploits  sparsity*  (and  the  variable  KEDGE  reduced  accordingly >,  or  other 
efficient  techniques  for  storing  sparse  tables  ma\'  be  used’'. 


C/ 

C/> 

C/> 

c/« 

c/> 


Subroutine  SETUP 

Purpose  :  SETUP  is  the  rr.ain  subroutine  to  generate  the 
eleinent-ele.’tient  adjacency  structure  from  the 


c/* 

element-node  incidence  table, 

INC,  for  a  planar 

* 

c/* 

finite  element  mesh. 

* 

c/« 

X 

0/* 

Input  Variables  : 

X 

c/» 

INC  -  element-node  incidence  table 

X 

0/  « 

X 

c/« 

Global  Variables  : 

X 

c/» 

NE  -  no.  of  elements  in  the  mesh 

X 

c/* 

NV  -  no.  of  nodes  in  the  mesh 

X 

c/« 

NVE  -  no.  of  nodes  oer  element 

X 

o/« 

NSPISI  -  (NE  -  1) 

X 

c/* 

NEDGE  -  possible  .maximum  no.  of  edge 

s  in  the  solution  graph 

X 

0/ * 

'(  =  NV  «  (NV-1)  /  2  ) 

X 

c/* 

EDGMAX  -  total  no.  of  edges  in  the  s 

elution  graph 

X 

0/ • 

ECGFLl  -  (EOGMAX  +  1) 

X 

c/» 

NETS  -  counter  for  the  edge-element 

structure  pair 

X 

c/» 

MAXS’JB  -  max.  subscript  used  in  the 

element -element 

X 

o/» 

adjacency  structure  pair  (see  subroutine  ELMAOJ) 

X 

c/» 

X 

c/» 

Output  : 

X 

0/ » 

The  element-element  adjacent  structu 

re  pair  (IA,CA)  is  stored 

X 

0/. 

in  the  array  PCCl  from  location  LOCI A 

to  (L0CJA*MAMSU5-1) . 

* 

0/. 

The  array  lA  is  stored  in  PCOLdCCIA) 

to  ?0CL(10CIA-NE) ; 

X 

c/» 

the  array  JA  is  stored  in  PCCL(LCCJA) 

to  ?CCL(LOCJA+MAXSUE-l) . 

X 

C/»  » 

C/»  Subroutines  used  :  * 

C/»  EDGSET,  EDG3LT,  ELKADJ  * 

C/» 

'2/ **1l******»****1l***1t*K1t1l*****1l*1t1tit1t1fK*1tnitnlfK1t*-KK**1f**1t1t**1C*9fK1t1C*1t**KnKit1( 

SUBROUTINE  SETUP (INC,  POOL) 

I.MPLICIT  INTEGER  (A-2) 

COKKON  /SYSTEM/  NE,  NV ,  NVE,  NEPLSl 
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C/- 

C/- 

C/- 

C/> 


c/« 

c/» 

c/« 

c/* 


COMMON  /GRAPH/  NEDGE,  EDGMAX,  EOGPLl ,  NELS,  MAXSUB 
COMMON  /lOSET/  IN,  10,  DEBUG 
EGG I CAL  DEBUG 

DIMENSION  INC(NE,NVE),  ?00L(1) 

Pnase  1  ;  zo  set  up  zhe  edge-element  adnacency  structure 
LOCNEL  =  1 

LOCSTA  =  LOCNEL  +  NEDGE 

CALL  EDGSETdNC,  POOL  (LOCNEL)  ,  POOL ( LOCSTA)  ) 

Phase  2  :  to  build  the  edge-element  adjacency 
structure  pair 

EDGPLl  =  EDGMAX  *  1 
LOCEL  =  LOCSTA  EDGPLl 
LCCPOS  =  LOCEL  -  NELS 

CALL  EDGBLTdNC,  POOL  { LOCNEL)  ,  POOL  (LOCSTA)  ,  POOL  (LOCEL), 
POCL(LOCPOS) ) 

Phase  3  :  to  build  the  element-element  adjacency 
structure  pair 


LOCI A  =  LCCPOS 
LOCJA  =  LCCIA  +  NEFL31 

CALL  ELKADJdNC,  POOL  (LOCNEL)  ,  POOL  (LOCSTA)  ,  ?CCL(LCCEL), 

*  POOL(LOCIA),  ?OOL(LCCJA)) 

RETURN 

END 

C/» 

C/»  Function  HASH 

C/*  Purpose  :  to  determine  the  location  of  an  edge  in 
C/»  the  symmetric  node  adjacency  matrix 

c/» 
c/* 

C/o 

c/« 
c/« 
c/* 
c/« 
c/* 

INTEGER  FUNCTION  HASH(N1,N2) 

IMPLICIT  INTEGER  (A-Z) 

J  =  AMAXC(N1,N2) 

I  =  AMIN0(N1,N2) 

HASH  =  (  (J-2) * (J-1) )/2  -  I 

RETURN 

END 

C/»  * 

C/*  Subroutine  EDGSET  « 

Purpose  :  to  set  up  counter  START  cf  the  edge-element  » 

adjacency  structure  pair  (START,  EL)  * 

{EL  will  be  determined  in  subroutine  EDGELT.)  ♦ 


Input  Variables  : 

Nl,  N2  -  nodes  in  the  mesh  (or  graph) 

Output  Variable  : 

HASH  -  location  of  edge  (Ni,N2)  in  the  symmetric  node 
adjacency  matrix 


c/» 

c/» 

c/* 

c/» 

c/« 


Inout  Variable  : 
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C/»  INC  -  elemer.t-node  incidence  tacie 

C/« 

C/»  Cucpui  Variables  : 

C/»  NEL(HASKtNl,N2) )  -  edge  no.  assigned  to  node  pair  (N1,N2) 

C/*  {  Only  edges  with  more  than  one  incident  element 

C/*  are  assigned  a  numher.  } 

C/*  STAST  -  starting  position  for  the  edge-element  adjacency 

C/»  structure  pair 

C/»  {  STAHT()t+l)  -  STABTOc)  =  no.  of  elements  adjacent 

C/x  to  edge  )c  } 

C/« 

C/*  Intermediate  Variable  : 

C/*  NEL(HASK(N1,N2) )  -  in  step  1,  the  array  is  a  ncde-ad;acency 

C/»  matrix  denoting  nu.mfcer  of  elements  into  edge  (N1,N2; 

c/* 

SU3RCUTINE  EDGSETdNC,  NEL,  ST.ART) 

IMPLICir  I.VTSGER  (A-Z) 

COK.MON  /SYSTEK/  NE,  NV,  NVE,  NEPLSl 

COM.MON  /GRAPH/  NEDGE,  SDGMA.X,  EDGPLl,  NELS,  HAXSU3 

DIMENSION  INC(NE,NVE),  NEL(NEDGE),  START(l) 

C/x 

C/x  Initialization 

C/x 

DC  5  I  =  1, NEDGE 
NEL (I)  =  C 
5  CONTINUE 

C/x 

C/x  Step  1  :  to  build  the  node  adjacency  matrix 

C/x 

DO  10  I  =  1,NE 
NN  =  NVE 

DO  10  J  =  1,  NN-1 
N1  =  INC(I,J) 

IF(Nl.EQ.C)  GOTO  10 
DC  10  K  =  J+1,NN 
K2  =  INC(I,K) 

IF(N1.EQ.N2  .OR.  N2.EQ.0)  GOTO  1C 
EDGE  =  HASH{N1,N2) 

NEL (EDGE)  =  NEL (EDGE)  +  1 
10  CONTINUE 

C/x 

C/x  step  2  : 

C/x  (1)  to  set  up  counter  START  for  the 

C/x  edge-element  adjacency  structure 

C/x  (2)  to  assign  a  number  to  the  edges  (with  two 

C/x  or  more  incident  elements  ) . 

C/x 

POSIT  =  0 
NELS  =  1 

DC  15  I  =  1,  .NEDGE 

IF(NEL{I)  .LE.  1)  GOTO  14 
POSIT  =  POSIT  +  1 
START (POSIT)  =  NELS 
NELS  =  NELS  +  NEL (I) 

NEL(I)  =  POSIT 


NEL(I)  = 


14 


0 


15  CONTINUE 

EDGMAX  =  POSIT 

START  (EDGMA»1)  =  NELS 

RETURN 

END 

C/.. . . . . . ................ 

0/« 

C/«  Subroutine  EDGBLT 

C/«  Purpose  :  to  build  the  array  ED  of  the  edge-element  adjacency 
0/*  structure  pair  (START,  EL) . 

0/» 

C/»  Input  Variables  : 

0/«  INC  -  element-node  incidence  table 

0/*  NEL(lt)  -  edge  no.  assigned  to  node  pair  (N1,N2) 

C/*  {  It  =  .H.ASH(N1,N2)  } 

c/*  START  -  Starting  position  for  the  edge-element  ad3acency 

C/*  structure  pair 

0/*  (  START(lt+l)  -  START(li)  =  no.  of  elements  adjacent 

0/*  to  edge  )t  } 

C/« 

C/«  Output  Variable  : 

C/*  ED  -  edge-element  adjacency  structure  pair 

C/»  {  The  set  of  elements  incident  to  edge  K.  is  stored 

0/*  in  ED(START(lt)  )  to  ED(START(k*l) )  .  } 

C/* 

C/»  working  Variable  : 

C/»  POSIT(k)  -  pointer  to  current  empty  position  in  EL  for  edge  k 

c/« 

c/». 

SUBROUTINE  EDGBLTdNC,  NEL,  START,  EL,  POSIT) 

IMPLICIT  INTEGER  (A-2) 

COMMON  /SYSTEM/  NE,  NV,  NVE,  NEPL31 
COMMON  /GRAPH/  NEDGE,  EDGKAX,  EDGPLl ,  NEDS,  KA,HSuB 
DIMENSION  INC(NE,NVE),  NEL(NEDGE),  START ( EDGPLl ) ,  EL(NELS), 
POSIKEDGMAX) 

Initialize  working  vector  POSIT 

DO  10  I  =  1,  EDGMAX 
PCSIT(I)  =  START(I) 

CONTINUE 

Fill  edge-element  adjacency  EL 

DO  20  I  =  1,NE 
NN  =  NVE 

DC  20  J  =  1,  NN-1 
N1  =  INC(I,J) 

IF(Nl.EQ.O)  GOTO  20 
DO  20  K  =  J+1,  NN 
N2  =  INC(I,K) 

IF(N2.EQ.0  .OR.  N1.EQ.N2)  GOTO  20 
EDGE  =  HASH(N1,N2) 

IF (NEL (EDGE)  .EQ.  0)  GOTO  20 
EDGACT  =  NED (EDGE) 

EL(POSIT(EDGACT))  =  I 
POSIT (EDGACT)  =  POSIT (EDGACT)  *  1 


Inode  Ni: 


Inode  N2I 
I  check  self  locpi 
I  location  in  NED  I 

ledge  number  I 
I  assign  I  to  EDI 
I  update  PCS  I  Tier.  I 


C/* 

C/» 

C/* 


C/« 

0/* 

C/» 


20  CONTINUE 
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RETURN 

END 

C/. 

C/»  Subroutine  ELMADJ 

C/«  Purpose  :  to  set  up  tfie  element-element  adjacency  structure 

C/«  pair  (IADJ,JADJ). 

c/« 

C/»  Input  Variables  : 

C/»  INC  -  element-node  incidence  table 

C/»  NELOc)  -  edge  no.  assigned  to  node  pair  (N1,N2) 

C/*  {  )t  =  HASH(N1,N2)  } 

C/«  START  -  Starting  position  for  the  edge-element  adjacency 

C/*  structure  pair 

C/»  {  START(li+l)  -  STARTOt)  =  no.  of  elements  adjacent 

C/*  to  edge  k  ) 

C/»  EL  -  edge-element  adjacency  structure  pair 

C/«  {  The  set  of  elements  incident  to  edge  k  is  stored 

C/»  in  EL(START(k))  to  EL  (START  ()c+l) )  .  } 

C/« 

C/*  Output  Variables  : 

C/*  (lADJ.JADJ)  -  element-element  adjacency  structure  pair 

c/* 

SUBROUTINE  ELMADJ (INC,  NEL,  START,  EL,  lADJ,  JA2J) 

IMPLICIT  INTEGER  (A-2) 

COMMON  /SYSTEM/  NE,  NV,  NVE,  NEPLSl 
COMMON  /GRAPH/  HEDGE,  EDGMA.X,  EDGPLl ,  NSLS,  MAXSUB 
COM-MON  /lOSET/  IN,  10,  DEBUG 
LOGICAL  DEBUG 

DIMENSION  I.VC(NE,NVE)  ,  NEL(NEDGE),  START  (EDGPLl) ,  EL(NELS), 
»  lADJ  (NEPLSl)  ,  fADJd) 

C/* 

C/«  Initialize  lADJ(l) 

c/* 

lADJ(l)  =  1 

C/« 

C/*  For  each  element,  find  ail  incident  elements  a.nd 

C/«  aueue  them  in  JADJ 


DO  50  1  =  1,NE 
ELMENT  =  I 
LI  =  lADJ(I) 

FILL  =  0 
NN  =  NVE 

DO  40  J  =  1,  NN-1 
N1  =  INC(I,J) 

IF(Nl.ES.O)  GOTO  40 
DO  4C  K  =  J+1,  NN 
N2  =  INC(I,K) 

IF(N1.E2.N2  .OR.  N2.EQ.0)  GOTO  40 
EDGE  =  HASH(N1,N2) 

IF(NEL(EDGE) .EQ.C)  GOTO  40 
EDGNUM  =  NEL (EDGE) 

BEGIN  =  START (EDGNUM) 

END  =  START (SDGNUM-1)  -  1 
DO  30  lEL  =  BEGIN,  END 
ELM  =  EL(IEL) 


;for  element  i: 


:ncde  Ni: 


Inode  N2: 

:checM  self  loop: 
'.location  in  NELl 

ledge  number: 

:for  all  elements 
:  incident  to 
:  edge  EDGNUM: 

:a  neigh,  element. 
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IF(ELM.EQ.ELMENT)  GOTO  30 

L2  =  LI  +  FILL 

IF(FILL.EQ.O)  GOTO  20 

Isame  as  element  I 

DO  10  JJ  =  LI,  L2-1 

checlc  if  there 

IF(JADJ(JJ) .EQ.ELM)  GOTO  30 

is  a  duplication 

10 

CONTINUE 

in  list  : 

20 

JADJ(L2)  =  ELM 

FILL  =  FILL  +  1 

queue  it  in  JADJ: 

30 

CONTINUE 

40 

CONTINUE 

IADJ(ELMENT+1)  =  lADJ(ELMENT)  +  FILL 

;  update  lADJ! 

50 

C/» 

CONTINUE 

MAXSUB  =  IADJ(NEPL31) 

IF ( .NOT. DEBUG)  RETURN 

C/« 

C/* 

Print  out  results 

WRITE (10, 1100) 

DO  60  I  =  1,NE 

LI  =  lADJ(I) 

L2  =  IADJ(I■^1)  -  1 

IF(L2.LT.L1)  GOTO  60 

WRITE (10, 1000)  I,  LI,  L2,  (JADJ(K),  K=L1,L2) 

SO 

CONTINUE 

RETURN 

1000 

FOR.MAT(T5,  14,  TIC,  16,  T18,  16,  (T30,10I6)) 

1100 

FORMAT (IHl,  T2,  'ELEMENT',  TIO,  ’IDAJ(I)',  TIS, 

IADJ(I+1) ' , 

« 

T30,  'JADJ  (NEIGHBORING  ELEMENTS)  :',/) 

END 
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